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Missing fluctuations in masked regions in the sky can be reconstructed from fluctuations in the surround- 
ing unmasked regions if they are sufficiently smooth. We propose to reconstruct such missing fluctuations by 
iteratively applying a spherical harmonic expansion to ffuctuations in the unmasked region. The accuracy of 
reconstruction depends on the mask geometries, the spectrum of underlying density fluctuations, and the num- 
ber of iterations. For Gaussian fluctuations with the Harrison-Zel'dovich spectrum, our method provides more 
accurate restoration than naive methods using the brute-forth matrix inversion or the singular value decomposi- 
tion. 
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I. INTRODUCTION 

After the first data release of the cosmic microwave back- 
ground (CMB) temperature fluctuations observed by the 
Wilkinson Microwave Anisotropy Probe (WMAP), anoma- 
lous signatures in the CMB on large angular scales, so called 
"large-angle anomalies" have been claimed by many authors 
[1-11] and also confirmed recently by the Planck mission 
[12]. So far, the origin has been veiled in mystery. They may 
be due to (a) difference in a priori significance and a posteriori 
significance [13-16] (b) incomplete subtraction of foreground 
emissions [17-19] (c) contribution from large scale structures 
via the integrated Sachs-Wolfe effect [12, 20-28], or kinetic 
Sunyaev-Zel'dovich effect [29], (d) possible systematics from 
instruments [30] (e) incomplete treatment of masking [31], or 
(f) extensions of inflationary models [32-39]. 

In real observation, we cannot observe the whole sky with 
a sufficiently high signal-to-noise (S/N) ratio. A conservative 
approach is to mask out the region (e.g., the Zone of Avoid- 
ance), where the S/N ratio is low and just ignore the data on 
it. 

For estimating the power spectrum of fluctuations in the 
sky, we can use deconvolution techniques [e.g. 40] assuming 
an isotropy prior on the power spectrum. However, for esti- 
mating the density field itself, we need to develop methods 
that can reconstruct the phases (if expressed in complex num- 
bers) as well as the amplitudes of missing fluctuations. 

In order to reconstruct missing data on the masked region, 
we need to find the inverse of the masking operator. However, 
in general, the masking matrix is singular thus not invertible. 
We need to assume a certain prior information about the un- 
derlying data such as isotropy and smoothness [31, 41^3] or 
their derivatives [44] in order to regularize the inverse oper- 
ator. Since the result depends on the choice of the prior, the 
robustness of each reconstruction method should be mutually 
checked. 

In this paper, we propose a new method for regularizing 



the inverse of masking operators using the iterative spherical 
harmonic expansion (IHE). As our aim is to probe the origin of 
the large-angle CMB anomalies, we ignore the effect of noise 
in the signal in the following, which is a good approximation 
for large-angle fluctuations. 

This paper is organized as follows. In section II, we de- 
scribe the formulation of our reconstruction method using the 
iterative spherical harmonic expansion and our simulations for 
random Gaussian fluctuations. In section III, we show the ac- 
curacy of the reconstruction, which depends on the pixeliza- 
tion of the sky, the underlying power spectrum, and the num- 
ber of iteration. And then we compare the accuracy of our 
method to that of the brute-force and the singular value de- 
composition (SVD) method. In section IV, we give our con- 
clusion. 



II. ITERATIVE HARMONIC EXPANSION 
A. Theory 

Suppose a density field 6{y) on a unit sphere where y rep- 
resents a unit vector that specifies the angular position in the 
sky. The observed field is given by 



5obs(r) = mr)'^(r), 



(1) 
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where = in the masked region and W - \ elsewhere. 
Then we can obtain the pseudo harmonic coefficient 5, by sim- 
ply integrating over the whole sky, 

a,- = J'^" 5{yW{y)V{y) (2) 

where W,j = J t/O 'WY jY* is the (/, j) component of the mode 
coupling matrix due to incomplete sky coverage, 7,(7) is a 
spherical harmonic, and c/Q is the surface element on a unit 
sphere. For simplicity, we use the single subscript index / = 
-vl-\-m.+ \ where I and m represent a multipole moment and 
an azimuthal number, respectively. In terms of true harmonic 
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(a) Original map with mask 



(c) Pseudo-Qf,,, reconstmction 



(e) Pseudo-flf,,, residual 




(b) Original smoothed map without mask 



(d) IHE reconstruction 



(f) IHE residual 



FIG. 1: Plots of an original and reconstructed maps, (a) Original simulated map with 4ut = 30 and a sky cut \b\ < 20°, which is used as an 
input, (b) Original simulated map smoothed with 4ut = 7 on a full sky (no mask), (c) Reconstructed map with = 7 using pseudo-a,. 
(d) Reconstructed map with f^^^ = 7 using IHE with N^c = 10 iterations, (e) Residual map after subtraction: (b)-(c). (f) Residual map after 
subtraction:(d)-(b). Color scale (-1.87 to 1.55) is the same for all the figures. 



coefficients to be reconstructed a^'', Eq. (2) can be written as 

oo 

a, = «r^'V- (3) 
The mode coupling matrix W is explicitly given as [40] 

f(2A + l)(2^2 + l)(2^3 + l) 



= ;^«;,-,(-l)" 



£\ il ^3 1/ i 
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(4) 



where w, - j dQ WY* and (:::) denotes the Wigner 3j- 
symbols. Note that the subscript ij depends only on (j and 
nij, i.e. ij — + + f^j + 1- 

As shown in Eq. (3), af"^'s can be obtained by inverting 
the mode coupling matrix W. However, in general, it can be 
singular and not invertible. In order to regularize the inver- 
sion, we propose the following iterative scheme. The iteration 
process starts from 



dn 6(7Wif)Y;{fi 



(5) 



which we call "the pseudo-fl,'s". For m > I, the m-th set of 
a,'s is constructed iteratively from two maps; the original ob- 
served map at the unmasked region and the map reconstructed 
from the inverse transform of the (m - l)-th fl,'s at the masked 
region, 

af > ^ JdQ. [6(f)W(f) + ~6^"'-'\f)S(y)] y;(j), (6) 



where 



(7) = 



-(m-1) 







Yiif) m>2, 
m — 1, 



(7) 



and S - I - W. In what follows, we assume that fluctua- 
tions whose angular scales are larger than that of the masked 
region do not significantly correlate with fluctuations smaller 
than the masked region. In that case, the summation in Eq. (3) 
and Eq. (7) can be truncated at a certain multipole /,nax as long 
as we concern large-angle fluctuations corresponding to mul- 
tipoles i « /max- In what follows, we sum up fl,'s up to the 
multipole /jmax, and we omit the summation symbol when no 
confusion arises. 

Substituting Eq. (7) into (6) recursively, we obtain the gen- 
eral formula of the m-th iterated harmonic coefficients as a 
series of S a. 



~{m) 

a) 



= a]'\6fj + Sij + Si.Skj + ■■■ + [S'"-%), (8) 



where dfj is the Kronecker delta, and S ij - ^fj^ Wij. The map 

between m-th iterated harmonic coefficient 5*'"' is equivalent 

to the product of fl*'''s and the Taylor expansion of the inverse 
of the mode coupling matrix W truncated at the (m - l)-th 
order. 



a. ' 



(9) 



Thus a*"'*'s approximately represent fluctuations in which the 
mode coupling matrix W is deconvolved. The number of iter- 
ation A^ite should be evaluated using Monte-Carlo simulations, 
which is discussed later in Sec. III. 
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We can regard the above process as a mapping from a ob- 
servable to the estimator, 

FiHE : 5/ -> a^'"*. (10) 

As we will see in Sec. Ill, the accuracy of the IHE method 
depends on the maximum multipole (^^^ to be reconstructed, 
the spectrum of the underlying field 6, the mask geometry, 
and the number of iteration A^ite- Thus we can write Fihe - 

FlHEi{m..,6,W,Nn,). 

Note that the computation time for estimating 5*'^"=' can be 
significantly reduced if we use Eq. (6) instead of Eq. (8). The 
reason is as follows. As the rank of the matrix S ij is of the 
order t^^^^, it costs A^pix^max ^1- (4)' matrix algebra of 
Eq. (8) costs ~ Mte Anax Computations. If one uses Eq. (6), it 
would be of the order A^iteA^pix^max' where A^ite is the number 
of iteration, and A^pix is the number of pixels which tile the 
sky, and /"max is the maximum multipole to be reconstructed. 
For example, for given W = 10, A^pix = 12 X 1024- pixels 
and 5 time iteration, the computation time will be reduced by 
a factor of ~ 25. 
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B. Simulation 

In order to assess the accuracy of reconstruction, we gen- 
erate 1000 random isotropic Gaussian density fields in the 
sky. We use the code synfast, which is publicly available as 
a package in Healpix [50] [45] for generating random Gaus- 
sian maps. First, we use the Harrison-Zel'dovich spectrum as 
the input power spectrum. It gives an angular power spectrum 
Cf oc /" where n = -2 on large angular scales, in the Ein- 
stein de-Sitter universe which correspond to the Sachs-Wolfe 
plateau of the CMB power spectrum. As a simple model of 
the zone of avoidance, we consider an azimuthally symmet- 
ric mask W{y) = at the Galactic latitude h smaller than 20 
degrees, i.e. \h\ < 20° and W - I otherwise. We use pixels 
that are sufficiently smaller than the size of the mask and the 
fluctuation scales to be reconstructed in order to reduce the 
errors due to the pixelization effect (see Sec. Ill A). We adopt 
the Healpix resolution A^side =512 (the total number of pixels 
on the whole sky is A^pix = 12 x A^^.^^ = 3145728). The in- 
put power spectrum should be truncated at sufficiently small 
scales 4ut compared to the scale to be reconstructed ^niax, i-e. 
^max ^ 4ut- The cut-off scale 4ut is inferred from the de- 
tector resolution scale in CMB observations or the shot noise 
dominant scales in galaxy surveys. In our simulated maps, we 
set 4ut - 30, which is sufficiently smaller than the claimed 
anomalous multipoles in the CMB. We find that the choice of 
icut is not sensitive to our result. 



III. RESULT 

In this section, we first explore the pixelization effect, 
which may affect the accuracy of reconstruction in Sec. Ill A 
and then we compare our reconstruction method with other 
inversion methods using the brute-force or the singular value 



FIG. 2: The accuracy of map reconstruction averaged over 1000 
realizations of Gaussian maps as a function of the iteration number 
for a mask \b\ < 20° with fniax = 7 for various pixel resolution N^a^. 
Both solid and dashed lines are obtained by using Eqs. (6) and (7). 
Dashed lines represent values for which approximated 5 ,/s for pixel- 
based masks are used. Solid lines show values for which the exact 
expression of Eq. (11) is used. Dotted lines indicate values obtained 
by using the direct inversion of the mode coupling matrix W, the 
second equation in Eq. (9). In the bottom panel, we show the relative 
differences of the accuracy in percentage. A'sidc = 1024 provides 
accuracy better than a few percent. A'itc ~ 10 with A'sidc = 512 can 
achieve sub-percent accuracy. 



decomposition (SVD) in Sec. IIIB. Finally, we discuss the 
condition of the original maps that can be reconstructed by 
using the IHE in Sec.IIIC. 

A. Pixelization effect 

We need to assess the accuracy of the mode coupling ma- 
trix W since it may be close to being singular. In other words, 
the difference between W{f) and Wif,) should be carefully 
checked where f and 7, are unpixelized and pixelized posi- 
tions on a unit sphere, respectively. In order to do so, we need 
to estimate uncertainties due to pixelization of the sky as it 
could be potential sources of errors. For an azimuthally sym- 
metric mask W = for \b\ < bo, we do not need pixelization, 
since Wij can be analytically calculated as 

Wij = 2;r 1 - Yi(0,O) Yj(0, 0) sm0de\, (II) 

\ J,r/2-bo j 

where 6 is the polar angle, and 6 - njl - b. We compare 
the eigenvalues of W obtained by using Eq. (4) and Eq. (11) 
to study the accuracy of the calculated matrix. The differ- 
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ence depends on the resolution of pixels A^side- We adopt 
A^side = 16, 128,512, 1024 and find that the differences in the 
eigenvalues are less than 0.01% for A^side ^ 512. In order to 
evaluate the accuracy of map reconstruction, we use the 
norm of the fractional difference between the reconstructed 
map and the original map, 



(a) C( cc £0 



y'Vpix e2 



mi) 



(12) 



where 5rec and ^tme describe the reconstructed and the original 
density fluctuations contributed from the (-th multipole mode. 



S(f;{) = ^ a{„,Y[i„(f), 



(13) 



respectively. An overall fractional difference is defined by 
substituting 6{y) for 5{f, €) in Eq. (12) where 6{y) is a sum of 
all the multipole contributions, i.e. 5{y) - YjC 5iy\ 0- 

In Fig. 2, we can see the difference between a pixelized 
mask and an unpixelized "smooth" mask given by Eq. (11). 
Using Healpix, we pixelize the mask with different resolution 
A^side where the total number of pixel on the sky is defined as 



pix 



12 X A^^ijg. Before looking into the pixelization effect. 



we first study the overall behaviour of reconstruction accuracy 
as a function of the number of iteration. As we increase the 
number of iteration, the accuracy of reconstruction gradually 
improves at N-n^ < 10 then it begins to degrade. Finally, the 
D~ accuracy converges to the one obtained by using the brute- 
force inversion (dotted lines in the Fig. 2 ). As we can see in 
the bottom panel of Fig. 2, the differences in between the 
pixelized and unpixelized Wij are significant for A^side = 16 
and 128, while they are less than 1% at the minimum for 
A^side > 512. We then conclude that the pixelization effect can 
be ignored when we use sufficiently fine pixels with A^side ^ 
512. 



B. Comparison with ' and SVD 

Using Eq. (9), we can compare our method to the direct 
brute-force inversion of the matrix W and the singular value 
decomposition[46]. 

If the mode coupling matrix W is invertible, we obtain the 
unique solution of the underlying density fluctuation within 
the mask. However, if the matrix W contains some eigenval- 
ues which are close to zero, where the matrix is close to being 
singular, the inversion causes big errors. In such a case, we 
can remove the singularity by replacing the small eigenvalues 
with zero, so called the SVD method. More specifically, the 
matrix W can be decomposed into three matrices as. 



(14) 



where V and U are /^ax x 'max unitary matrices and a super- 
script t denotes the Hermitian conjugate. S is the diagonal 
matrix which consists of the eigenvalues of W, or called as the 
singular values. The order of the eigenvalues in 2 is arbitrary 
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FIG. 3; In the right panels, we plot accuracy as a function of 
iteration number. Three different rows show the different underlying 
density fluctuations that have power-law spectrum of Ct oz f, {^^ 
and from top to bottom. For a red (n = 0) spectrum, we do 
not see any improvement. However, for a flat (n = -2) spectrum, 
there exists an optimal number of iteration around A^i^ = 10 that 
minimizes D". The optimal N^c depends on the For a blue 

(n = -4) spectrum, D- converges to a few percent accuracy after 
sufficient number of iteration. For all cases, after a sufficient number 
of iteration, converges to the value which is obtained by the direct 
inversion. In the left panels, we show the result for the SVD method 
for the comparison as a function of eigenvalue threshold Aj,. 



but they are arranged in a descending order so that the decom- 
position is determined uniquely. Let the A:-th eigenvalue be 
Af,, and the eigenvalues smaller than A/^ are set to zero. Then 
the pseudo inversion of W is written in terms of the rank-^ 
diagonal matrix that consists of the reciprocal of the non-zero 
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FIG. 4: The IHE reconstruction accuracy D-^ for each multipole component. From right to left, the indices of underlying power spectrum are 
n = 0, -2 and -4, respectively. From top to bottom, the components of multipoles = 1 to 4 are shown. We reconstruct the masked map with 
maximum multipoles of ^max = 3, 5, 7 and 9. The error bars represent 1 cr standard deviation calculated from 1000 random realizations. 



eigenvalues that are equal to or larger than A^, as 

W*^VL^U\ (15) 

which gives 

af = hjWtj. (16) 

The choice of threshold is not trivial, and should be care- 
fully determined a priori because the mapping Fsvd depends 
on various factors including the threshold Ak'. i.e. Fsvd = 
■P'svD(^max, <5, W, Ak). In this context, the singularity can only 
be defined in terms of A^. More specifically, the mapping 
FsvD is contaminated by the eigenvalues which are close to 
being singular if dD'/dAt < 0. In such cases, the threshold 
of eigenvalues should be increased to implement more accu- 
rate reconstruction. The brute-force inversion corresponds to 

^'inv = ^'sVD( ^max, S, W, Ak ^ 0). 

In the left panels of Fig. 3, D~ is shown as a function of 
the minimum non-zero eigenvalue A^. The solid lines show 
the different maximum multipole to be reconstructed, f^^^ - 
3,5,7,9 from dark to light color. The different three rows 
show the different power of the input power spectra which 
will be discussed in detail in the next section. Let us focus on 
the middle panel. For ^^ax = 3 and 5, D~ increases monoton- 
ically as we increase the threshold. It means that the Fsvd is 
not contaminated by singular eigenvalues as dD'/dA^ > 0. On 
the other hand, for ^^ax = 7 and 9, has a minimum around 
Aii ~ OA. Thus the eigenvalues smaller than ~ 0. 1 may be the 
contaminant of the mapping F and we can better estimate the 
original density fluctuations when we limit the eigenvalues as 
A < 0.1. In the right panel of Fig. 3, we show D~ for our IHE 
method as a function of the number of iteration. As in Fig. 2, 
has the minimum around A^ite - 10 depending on the max- 
imum multipole. After a sufficient number of iteration, D~ 
values converge to those obtained by using the brute-force in- 
version, which are shown as the horizontal dashed lines in 



Fig. 2. Our IHE method with a certain finite number of itera- 
tion is always better than the SVD method. However, in prac- 
tice, we should know a priori the optimal number of iteration. 
The number depends on the mask geometry, the maximum 
multipole to be reconstructed, and the underlying spectrum of 
the density fluctuation. We can estimate the optimal iteration 
number by carrying out Monte-Carlo simulations. In the next 
section we will see how the result changes for different types 
of underlying power spectrum. 

C. Dependence on underlying power spectrum 

Reconstruction accuracy depends on the underlying power 
spectrum Cf oc i" of the fluctuation as each harmonic mode 
is not independent in the masked incomplete sky even the un- 
derlying fluctuation is Gaussian. We consider three power law 
indices, n - 0, -2 and -4, which correspond to the following 
three cases. The projected two dimensional galaxy or dark 
matter distribution is approximated as oc [e.g. 47] on 
large scales, and the ordinary Sachs-Wolfe spectrum is ap- 
proximated as C^* oc [48]. On very large-angular scales, 
the integrated Sachs-Wolfe effect, which gives CJ?* oc (^'^ 
[e.g. 49] dominates the CMB in the standard ACDM scenario. 

Given the typical scale of mask 6?m, it is impossible to re- 
construct the fluctuation whose scale is smaller than (y[ > 
18O/0M- Because of the mode coupling, fluctuations with an- 
gular scales corresponding to /'m are strongly affected by fluc- 
tuations with smaller angular sizes. If the spectral index is 
negative, the amplitude of a smaller scale fluctuation is weak 
and does not strongly disturb the large scale fluctuations. Thus 
the deconvolution mapping F is less affected by singularities. 
However, if the slope of the spectrum is flat or positive, large 
scale modes are highly contaminated by the noisy small scale 
fluctuations. 

In Fig. 3, the top and bottom panels show accuracy for 
« = and « - -4. For « = case, the SVD has a min- 
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imum at Ak ~ 0.5 that is larger than the case of « = -2. 
It means that Fsvd(« - 0) is more affected by singularities 
than FsY£,(n - -2). The IHE result in the right panels shows 
A^ite = 1 gives the best accuracy and it gradually degrades to 
the value given by the brute-force method. On the other hand, 
Fsvoin - -4) is not affected by singularities because it al- 
ways shows dD^/dAk < for the SVD method. In Fig. 4, 
we also show the reconstruction accuracy for each multipole 
component, Dj. The behaviour is similar to the one described 
in Fig. 3. The reconstruction accuracy is always better for low 
multipoles and that the optimal number of iteration depends 
also on the multipole component we are going to reconstruct. 

IV. SUMMARY 

We have developed a new method based on the iterative har- 
monic expansions (IHE) for reconstructing the missing fluc- 
tuations in the masked region in the sky. Reconstructing the 
data at a masked region is known as an inverse problem. Our 
method is equivalent to the brute-force inversion in the limit 
that the iteration number goes to infinity. However, in some 
cases, the finite truncation of the iteration gives a better esti- 
mate of the underlying fluctuation. The reconstruction accu- 
racy depends on the geometrical shape of the mask, the max- 
imum multipole mode to be included in the IHE analysis, the 
tilt of the underlying power spectrum, and the multipole com- 
ponent of the map to be reconstructed. 

As an example, we have applied the method to reconstruct 



the missing data on an azimuthally symmetric mask. We 
have considered three types of gaussian fluctuations with the 
power-law indices n - 0, -2 and -4, which correspond to the 
galaxy power spectrum, the ordinary Sachs-Wolfe spectrum, 
and the integrated Sachs-Wolfe spectrum, respectively. For 
« = -2 case, we have found that there exists an optimal fi- 
nite number of iteration that makes the reconstruction more 
accurate than the SVD method or the brute-force method. For 
« = case, the pseudo-af,,, is the best estimator for the pro- 
jected density fluctuations. For n - —A case, the brute-force 
inversion method gives the best accuracy. In that case, the IHE 
method can help to reduce the computation time of inversion. 

It would be interesting to see how the significance of the 
large-angle CMB anomaly changes when we use different 
methods of map reconstruction. The IHE method may help 
to unveil the origin of the CMB anomaly. We will explore this 
problem in our future work. 
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